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Abstract 


In Irrigation District No. 38, Río Mayo, in the state of Sonora, México, 
wheat cultivation represents around 70 % of the area under irrigation 
(around 68,000 ha), so several studies have been performed with the 
aim of relating the crop yields with values of the NDVI index obtained 
from images by satellites Landsat 7, Landsat 8 and Deimos, as well as 
with estimations of evapotranspiration of the crops (ETc), evaluated 
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through two methods, with the ETc = K¿p*ETr relation (ETr, Reference 
evapotranspiration), where K., is estimated as a function of the NVDI 
index, and through the MATRIC moisture balance method. Specifically, 
for the early estimation, a relationship between the yield and the value 
obtained with the NDVI index during the crop's flowering has been 
considered, a time when it generally achieves its highest value and 
because there are enough antecedents in the literature, among them 
studies performed by the International Water Management Institute and 
“Colegio de Postgraduados”, which stand out. However, this NDVI value 
is not enough for the estimation, so the value accumulated from heat 
units can be used, as will be shown in this article. On the other hand, in 
order to communicate the conditions of their crops to agricultural 
producers, viewers have been built, accessible through the Internet or 
local networks. Also, as part of the PLEIADeS and SIRIUS projects, 
financed by the European Commission, where the authors collaborate, 
the SPIDER viewer was created, which allows observing the NDVI and 
ETc values in each plot, although there are others, one designed by one 
of the authors and another by a doctoral student from Colegio de 
Postgraduados; the latter is installed in the Irrigation District. To 
calibrate the methodology proposed in estimating the probable yield, 
operations staff from the irrigation district performed measurements of 
the yields in several plots, with the objective of finding a function 
between the yields and the NDVI values during the crop's flowering. 
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Background 


In Report No. 28 of the International Water Management Institute 
(IWMI), related with the evaluation of irrigated area in the Bhakra 
irrigation system in India (R. Sakthivadivel et a/., 1999), using remote 
sensors and geographic information systems (GIS), a classification of 
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the productivity of the land is made by using satellite images with which 
the values of the Normalized Difference Vegetation Index (NDVI) are 
calculated, and they are related with yields from the wheat crop, which 
dominates this irrigation system of 1.3 million hectares. 


The basis of this evaluation is the relationship between the fraction of 
the photosynthetically active radiation and the crop's yield. Indeed, this 
radiation called PAR describes the radiation available for the vegetation's 
photosynthesis process. On the other hand, the fraction of solar 
radiation absorbed by the chlorophyll pigments (fPAR) describes the 
energy related with the assimilation of carbon dioxide and is derived 
from the PAR absorbed by the vegetative canopy, divided by the PAR 
available from sun radiation. It is assumed that fPAR is worth zero for 
the bare soil. 


In this report, it is assumed that the photosynthetically active radiation 
absorbed by the crop (APAR) is related with the sum in time of the fPAR 
whose units are: (J] m? t?). It is also pointed out that APAR is the main 
parameter that controls the total biomass accumulated by the crop, so 
that the yield can be present as the following function: 


Y =zNAPAR (1) 


Where Y is the grain yield in the crop and z is the relationship between 
the grain yield and the total biomass, and »7(kg/J) is the efficiency of 
photosynthesis depending on the type of route for carbon fixation (C?, 
C*, or CAM?) (Prince 1991). The validity of the equation (1) has been 
proven by Monteith (1972), and experimentally by Daughtry et al. 
(1992), Field, Randerson and Malmstorm (1995) and Morel et al (2014). 


Then: —fPAR =a+b*NDVI, (2) 


Where a and b are constants and NDVI is an expression directly linked 
with the crop's yield. Based on this report, there is reference to several 
studies where the value of NDVI in anthesis (crop heading stage) is 
related, in order to estimate the crop's yield, because in this stage an 
evaluation of the final yield of the crop can be made; in this regard, they 
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mention the study by Pestamalci et al. (1995). However, in México a 
good linear relation was obtained between the yield in barley and the 
NDVI value in this stage of the crop's development (Ruiz Huanca Paulino 
et al., 2005). 


For the case of wheat in India, the function found was the following: 
Y = 10.99 NDVI - 3.75 (3) 


With R* = 0.8594 and standard error of 0.217 t/ha. 
For the case of barley in México, the equation found was: 


Y = 11.48 NDVI -0.882 (4) 


with R? = 0.8846 


In México also, a function for wheat in the south of the state of Sonora 
has been found, with a similar structure: 


Y = 12.88 NDVI - 4.13 (5) 


with R? = 0.873 and standard error of 0.495 t/ha. 


Materials and methods 


The study research was performed in Irrigation District number 038, 
“Río Mayo”, in the south of the state of Sonora, located at coordinates - 
109.7873 W, 27.2912 N upper left corner and -109.2159 W and 
26.7059 N lower right corner; it is irrigated with water from the "Adolfo 
Ruiz Cortínez” dam and 130 wells. The irrigable surface is 98,520 ha. 


In this district, the geographic information system (GIS) is available 
with 15,022 plots detected, and it is divided into 16 irrigation modules, 
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which are managed by 16 Water Users” Associations (WUA) and a 
Limited Liability Society, which oversees the operation and maintenance 
of the infrastructure of the major network of canals, drainage and 
structures. 


During the last four agricultural cycles, follow-up of the surface irrigated 
has been performed through monitoring with satellite images by Landsat 
7, Landsat 8, and in a few months with Deimos 1 and RapidEye. To 
perform this monitoring, the NDVI (normalized difference vegetation 
index, Rouse et al., 1974) has been calculated, and the MSI (moisture 
stress index, Rock et al., 1986), using the corresponding reflectance 
calculated from satellite images, to which an atmospheric correction was 
applied and a  multi-sensor  normalization. In addition, the 
evapotranspiration of crops has been estimated as: 


ETc = Kc*ETr (6) 


Where ETc is the crop's evaporation and ETr is the evapotranspiration of 
reference, generally estimated through the Penman-Monteith equation, 
and the calculation of basal Kc as an empirical linear function of the 
NDVI (D'Urso and Calera, 2005, Calera et al, 2017). 


The functions mentioned are: 


NDVI = Pi Pr (7) 
PFP, 


Where: p; close infrared reflectance and p, reflectance in red 


mMsI =P (8) 
Pim 


Where: p; close infrared reflectance and Pjm mid-infrared. 


Kc = 1.26 NDVI - 0.2 (9) 
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Coefficients of this function can vary depending on the crop; the values 
shown are valid for the wheat crop. 
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crop's ETc for each 


plot. 


The average values of the indexes for each plot are estimated using a 
program developed by one of the authors of this article, and it calculates 
the average value of pixels plus fraction, within the limits of each GIS 
plot. 


For the evaluation at the level of module or plot, modeling of the NDVI 
development (accumulated values) was carried out, assuming that this 
growth is proportional to the growth of the biomass and that it adjusts 
well to a sigmoid model, as J. H. Thornley (1976) suggests. Generally, 
between these models, those most frequently used are the Logistic and 
the Gompertz (Yin et al, 2003). 


The differential equations on which the biomass development is based, 
for the case of the Logistic model, is 
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d 
LS = E, —y)y (10) 
dt 


For the Gompertz model it ¡s: 


S ==k|LnY, — Lny)y (11) 


Where y is the dependent variable which, for this case, could be the 
amount of biomass, t the development time, k a constant of 
proportionality, and Y,, the maximum value that the y variable reaches. 


When integrating, the following functions are obtained: 


y = — 
14 Ln 0 exp(—ct) 


Yo 


Which can appear as 


a 
s 12 
e 1+bexp(-ct) (12) 


Yo is the initial value of the variable, which from (12) for t = O is worth: 


a a 
1+b 


Yo 
For the Gompertz function it ¡s: 


y = Y, exp[-exp(b—ke)] 


Which usually appears as: 
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$e aexp|[-exp(b—ct)] (13) 


The growth rate is also usually used in both cases, based on functions 
(10) and (11); these rates that are presented as a bell function have a 
maximum value at the inflection point and are a good indicator of when 
the maturation and senescence of the crop begin. For the case of the 
Logistic function, the maximum value is reached at a time 


17 (14) 


For the case of the Gompertz function, this value is reached when 


el 115) 


m 


Specifically, for the case in consideration, it was seen that a better 
adjustment was achieved with the Logistic function, usually with values 
of the determination coefficient R? higher than 0.99 and standard errors 
below 0.205 in the NDVI. 


With regard to access to graphic satellite information by users of the Río 
Mayo irrigation district, they can enter the Web and see the state of 
their plots. First, they could access the information through the viewer 
generated in the PLEIADeS and SIRIUS projects, with servers in Spain, 
but since these projects are finished, another one was developed by 
“Colegio de Postgraduados” as part of a Master's thesis, which is 
installed in the irrigation district of study. 


Finally, in order to have another tool that allows performing an early 
estimation of yield, the heat units (Growing Degree Days: GDD) in 
Celsius degrees will be used, which allow evaluating the energy 
accumulated, from sowing to flowering, as well as the total accumulated 
until the end of maturation. 
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The climate information is obtained from INIFAP's automatic station in 
the irrigation district's center, for the whole period of crop development, 
daily, from October to May, by the state network of meteorological 
stations (Patronato para la Investigación y Experimentación Agrícola del 
Estado de Sonora) at: http://pieaes.dyndns.org/. 


Results 


During the last four agricultural cycles: 2011-2012, 2012-2013, 2013- 
2014 and 2014-2015, measurements were made of the wheat yield 
through sampling in several of the district's irrigation modules, where 
this crop predominates, occupying between 75 % and 80 % of the 
module's surface. In the first two cycles, a good relation was found 
between the maximum level of the NDVI and the average wheat yield; 
however, in the 2013-2014 cycle a small decrease was observed in the 
yield with regard to the expected one in agreement with function (5), 
but for the 2014-2015 agricultural cycle, the decrease in yield was 
highly significant, and the differences estimated in 6 of the modules 
evaluated varied between 20 % and 30 % less than in 2013-2014. 


For the case of module 15, the estimation of the average wheat yield in 
the module was obtained through average value sampling of 6.91 t/ha 
in 2013, 6.62 t/ha in 2014 and 5.19 t/ha in 2015, a reduction of 1.72 
t/ha equivalent almost to 25% of the yield obtained in 2013. 


To understand the differences observed, the average NDVI values 
observed were graphed against the Heat Units: Growing Degree Days 
(9C) and the values for the maximum NDVI each year were also 
graphed against the yields as shown in Figure 2 (Parthasarathi et al, 
2013). 
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Figure 2. Comparison of the NDVI values and yields against heat units. 


As complement to this figure, Table 1 shows a summary of the 
information that can be obtained with it. 


Table 1. Information on Module 15 


Max GDD. Yield. Max Day GDD Max 
Year NDVI total t/ha NDVI Val. 
2013 0.89 2331 6.91 128 1663 
2014 0.73 2169 6.62 144 1573 
2015 0.77 1809 5.19 96 1180 


From Figure 2 and Table 1, it is deduced that the accumulation of Heat 
Units when the maximum value of the NDVI was present was reduced 
significantly in 2015 as compared to the prior years, and the same 
happens with the total HUs accumulated at the end of the crop, when 
the NDVI values reach a minimum of 0.3. These decreases imply that 
the crop did not accumulate enough energy to attain higher yields. 


An analysis of the crop's behavior was also performed in other modules 
where wheat is the predominant crop, Modules 5, 6, 10 and 13, whose 
results are shown summarized in Table 2 presented next. 
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Table 2. Information about the yield variation in 5 Modules. 
GDD at 
NDVI Total 
Module Year max GDD Yield t/ha Yield dif. 
Mod_05 2014 1573 2169 6.24 
Mod_05 2015 1006 1678 4.81 1.43 
Mod_06 2014 1573 2463 6.51 
Mod_06 2015 1180 1809 5,42 1.09 
Mod_10 2014 2046 2464 7.08 
Mod_10 2015 1180 1921 5.19 1.89 
Mod_13 2014 1573 2464 6.83 
Mod_13 2015 1180 1809 4.68 2.15 
Mod_15 2014 1573 2169 6.62 
Mod_15 2015 1180 1809 5.19 1.43 


As can be observed in this table, there was a significant reduction of the 
yields in every case, in all the irrigation Modules considered; in the case 
of Module 13, the reduction was almost a third of the yield obtained in 
2014. 


In general, there was a significant reduction of the Heat Units, which 
indicates a lower accumulation of energy in the crop, generating the 
yield reduction. To analyze the yield reduction in more detail, a more 
specific analysis for the case of Module 15 has been performed, where 
there is good information about the yields because the sampling was 
broader and more careful as a result of the participation of its directors 
who are producers that rent considerable surfaces. 


For this case, modeling of the crop growth was performed, using the 
NDVI indexes accumulated, adjusting the variation to a sigmoid function 
of the Logistic type, and the rates of wheat crop growth were calculated 
for each one of the three years when their yields were available. 


Adjustment of the values accumulated to a Logistic or Gompertz 
function was tested. The data were adjusted better to the Logistic 
model. 
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For 2013, the model: y= Esad with R? = 0.9916 
1+33.638 * exp(-0.0298 * +1) 


Standard Error= 0.202 


For 2014, the model: y = E 0Te with R? = 0.994 
1+40.137 *exp(-0.035*1) 


Standard Error =0.128 


For 2015, the model: y= 05 with R? = 0.989 
1+11.769*exp(-0.032*1) 


Standard Error =0.201 


In Figure 3, the curves generated with these functions have been 
graphed, showing that the longer the crop lasts, there was higher yield; 
it should be noted that in 2013, there was a longer duration of the crop 
than in 2015, and in 2014 it can be appreciated that there is an 
intermediate value. 


Growth models in 3 years 
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Figure 3. Adjustment of a Logistic function of accumulated NDVI 
values. 


When calculating the growth rates of the wheat crop in these three 
years, it is observed that the values of the Heat Units for maximum 
NDVI value each year have differences, which are accentuated 
significantly between 2013 and 2015, as is shown in Figure 4, where in 
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addition to graphing the rates with regard to the Hus, it has also been 
done in days. In both cases, the difference between 2013 and 2015 can 
be noted. 


> > 
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Figure 4. Growth rates of the NDVI accumulated in three years, Module 
15. 


In the figures it can be seen that the growth rates in years 2013 and 
2014 show certain similarities, with the variant of a higher maximum 
value in 2013; however, it is evident that the 2015 curve is completely 
out of phase, the maximum value happened prematurely, which agrees 
with the position of the NDVI values of the images shown in Figure 2. 


The reason of the notable variation observed in the year 2015, among 
other variant meteorological conditions, was observed in this year, with 
many cloudy days, which probably generated a greenhouse effect, 
accelerating the development of the crop and preventing energy 
accumulation from happening, which is noted in the lower accumulation 
of HU. 


When calculating the average evapotranspiration for module 15, it is 
observed that in spite of the differences in the crop's development, the 
evapotranspiration in year 2013 was virtually equal to that in 2015, of 
around 390 mm, which implies that the daily evapotranspiration rate in 
2015 was higher; this fact can also explain the greenhouse effect 
generated by cloudiness, which accelerated the crop's growth. 
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It can be concluded that under “normal” meteorological conditions, an 
estimate of the yield value can be made based on the maximum NDVI 
value observed, which generally agrees with the time when anthesis 
happens, or flowering of the crop, as has been shown in the studies 
developed in different parts of the world, among which it can been seen 
in India and México. 


However, a very fast growth of the biomass does not allow accumulating 
enough energy and, therefore, a decrease in yield can be expected, 
which could be proportional to the deviation observed in HU and even 
simply in days, with regard to the more frequent or “normal” value. 
Thus, an early estimation of the yield could be performed based on the 
maximum level of the NDVI, corrected in function of the date of its 
presentation, which allows evaluating the heat units accumulated. 


Naturally, this work will require testing in other sites and possibly with 
more measurements to have a greater certainty about these 
conclusions. It is worth mentioning that there is more detailed 
information regarding the yield values in several plots, primarily in 
irrigation modules 13 and 15, with which a thesis study is being carried 
out to corroborate the prior results that are presented in this article. 
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